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ABSTRACT 

Self-gravitating bosonic fields can support stable and localised (solitonic) field 
configurations. Such solitons should be ubiquitous in models of axion dark 
matter, with their characteristic mass and size depending on some inverse 
power of the axion mass, rua- Using a scaling symmetry and the uncertainty 
principle, the soliton core size can be related to the central density and ax¬ 
ion mass in a universal way. Solitons have a constant central density due to 
pressure-support, unlike the cuspy profile of cold dark matter (CDM). Conse¬ 
quently, solitons composed of ultra-light axions (ULAs) may resolve the ‘cusp- 
core’ problem of CDM. In DM halos, thermodynamics will lead to a CDM-like 
Navarro-Frenk-White (NFW) profile at large radii, with a central soliton core 
at small radii. Using Monte-Carlo techniques to explore the possible density 
profiles of this form, a fit to stellar-kinematical data of dwarf spheroidal galax¬ 
ies is performed. The data favour cores, and show no preference concerning 
the NFW part of the halo. In order for ULAs to resolve the cusp-core problem 
(without recourse to baryon feedback, or other astrophysical effects) the axion 
mass must satisfy rUa < 1.1 x 10“^^ eV at 95% C.L. However, ULAs with 
1 X 10“^^ eV are in some tension with cosmological structure forma¬ 
tion. An axion solution to the cusp-core problem thus makes novel predictions 
for future measurements of the epoch of reionisation. On the other hand, im¬ 
proved measurements of structure formation could soon impose a Catch 22 on 
axion/scalar field DM, similar to the case of warm DM. 

Key words: Cosmology: theory, dark matter, elementary particles - galaxies: 
dwarf, halos. 


1 INTRODUCTION 


Dark matter (DM) is known to comprise the majority of 


the matter content of the universe (e.g. Planck Collabo- 
|ration|2014[ |2015[ ) . The simplest and leading candidate 
is cold (C)DM. CDM has vanishing equation of state 
and sound speed, w = Cs = 0, and clusters on all scales. 
Popular CDM candidates are ©(GeV) mass thermally 
produced supersymmetric weakly interacting massive 
particles (SUSY WIMPs, e.g. Jungman et al.|19^ l, and 
the 0{fieY) mass non-thermally produced QCD axion 


(Peccei & Quinn 1977 Weinberg]|1978 Wilczek 19781. 
The free-streaming and decoupling lengths of a WIMP, 
and the Jeans scale of the QCD axion are both ex- 
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tremely small (i.e. sub-solar on a mass scale, see e.g. 
Loeb fc Zaldarriaga||2005 f 


It is well known, however, that CDM faces a num¬ 
ber of ‘small-scale’ problems related to galaxy forma¬ 
tion: ‘missing satellites’ | Moore et al.[199'9 Klypin et al. 


19991, ‘too-big-to-fail’ ( Boylan-Kolchin et al.|20TT I, and 


cusp-core’ (reviewed in|Wyse fc Gilmore 20081, to name 


a few. Baryonic physics, for example feedback or dynam¬ 
ical friction, offers possibilities to resolve some or all of 


these problems within the framework of CDM (e.g. 

Del 

Popolo 

120091 

Governato et al. 20121 Pontzen & Gov- 

ernatol 

2014 1 

3el Popolo & Pace 2015| and references 


therein). Modifying the particle physics of DM so that 
it is no longer cold and collisionless also provides an at¬ 


tractive and competitive solution (e.g. Spergel & Stein- 


hardt 

|20001 Hu et al. 1120001 Peebles 

to 

o 

o 

o 

Boehm et al. 

200 1| 

Bode et al.|200111Marsh & Sil 

k|2013 

lElbert et al. 
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20141. Modified gravity has also been considered in this 
context (e.g. Lombriser &: Penarrubia]|2014 1. 

There is no compelling theoretical reason, however, 
that the DM should be cold or indeed thermal. From a 
pragmatic point of view, then, possible signatures of the 
particle nature of DM on galactic scales provide a use¬ 
ful tool to constrain or exclude models. The small-scale 
problems provide a useful way to frame our tests of DM 
using clustering. If the small-scale problems are resolved 
by baryon physics, then this removes one ‘side’ of par¬ 
ticle physics constraints derived from them. However, 
the other side is left intact: a particle physics solution 
cannot ‘over solve’ the problem (e.g. too little structure 
or too large cores). 

Thermal velocities and degeneracy pressure can 
support cores if the DM is warm (WDM, e.g. [Bond et al. 
1982 Bode et al.||2001 l. However, mass ranges allowed 
by large-scale structure constraints, mw ^ ©(few keV) 


limit core sizes to be too small (e.g. Maccio et al. 


2012 Schneider et al. 20141. Particle physics candi¬ 


dates for WDM include sterile neutrinos, or the grav- 
itino. Another promising solution is offered by ul tra- 
light (pseudo) scalar field DM with rUa ~ 10“^^ eV (Hu 
et al.|[?000[ [Marsh fc Silk 2013 Schive et al. 2014a I 
which may be in the form of axions arising in string 
theory ( Svrcek fc WittelilpOOO Arvanitaki et al)||2010 ) 
or other extensions of the standard model of particle 


physics (see Kim 1987 for a historic review). In this 
model, cores are supported by quantum pressure aris¬ 
ing from the axion de Broglie wavelength. Axions with 
TUa > 10“^® eV are consistent with large-scale structure 
constraints (Bozek et al. 2014|) and can provide large 


cores ijHu et al.||2000||Marsh fc Silk||2013| [Schive et al. 


2014a|b I. Such ultra-light axions (ULAs) may therefore 

offer a viable particle physics solution to the small-scale 
problems. 

If the small-scale problems are in fact resolved by 
ULAs, then forthcoming experiments, such as ACTPol 


(Calabrese et al. 20141, the James Webb Space Tele- 
scope (JW ST, [Windhorst et al.|2006p, a nd Euclid ( [Lau-| 
reijs et n)|[2011[ ' Amendola et al.|[2013 l, will find novel 
signatures incompatible with CDM. These include a 
truncated and delayed reionisation history ([Bozek et al.| 
20141, a dearth of high-a galaxies (Bozek et al. 20141, 


and a lack of weak-lensing shear-power on small scales 
(by analogy to WDM, e.g. [Smith fc Markovic[[20TT| . If 
these observations are consistent with CDM, however, 
then ULAs will be excluded from any relevant role in 
the small-scale crises. 

We note that there is not consensus on the prefer¬ 
ence for cores versus cusps in dwarf density profiles (e.g. 


Breddels fc Helnn|[2013[ [Richardson fc Fairbairn||2014[ 


Strigari et al 20141. If dwarfs are in fact cuspy, then 


no baryon feedback or dynamical friction is necessary, 
and light DM providing cores would be excluded. In this 


work, we use the data of Walker & Penarrubia (20111 


(hereafter, WPll), who report that their measurements 
exclude NFW-like cuspy profiles at a confidence level of 
95.9% for Fornax and 99.8% for Sculptor. 


The nature of axion clustering on small scales is 
a fascinating topic regardless of their role or otherwise 
in the small-scale crises. There has been much discus¬ 
sion in the literature (e.g. [Sikivie fc Yang[[2009 


David- 


son fc Elmer[[2013[ [Davidson[[2015[ [Guth et al.[[2014| 


Banik fc Sikivie[[2015[ and references therein) concern¬ 


ing whether axions, including the QCD axion, undergo 
Bose-Einstein condensation and display long-range cor¬ 
relation. This question is of more than theoretical im¬ 
portance and can greatly affect the direct detection 
prospects for the axion, for example by ADMX puffyl 


et al. 2006 Hoskins et al. 20111. The density solitons 


that we study here are the same ground-state solutions 
studied in the context of the QCD axion by [Guth et al.[ 
(20141, and display only short-range order. For ULAs, 


these solitons are kpc scale, while for the QCD axion 
they are closer to km scale. 

DM composed of axions of any mass, be it the 
QCD axion or ULAs, possesses a characteristic scale and 


DM clustering should be granular on some scale (Schive 
et al.|2014al. This scale is set by the axion Jeans scale. 


which varies with cosmic time, making structure forma¬ 
tion non-hierarchical (Marsh fc Silk[[2013 Bozek et al. 


20141. This departure from CDM on small scales makes 


the study of axion clustering on these scales impossi¬ 
ble using standard N-body techniques, and we must ask 
many basic questions afresh. What fraction of the DM 
in the Milky Way is smooth, and what is in solitons? 
What is the power spectrum on small scales? What is 
the mass function of solitons and sub-halos? Answering 
these questions in detail will be the subject of future 
work, with the present work taking some small steps 
in this direction. The scale symmetry of the relevant 
equations for axion DM makes answers to these ques¬ 
tions universal, and equally applicable to ULAs and the 
QCD axion. 

In this paper, we study the soliton solutions of ax¬ 
ion DM, clarifying the core formation mechanism. The 
stability of this profile and its validity outside the spher¬ 
ically symmetric case is supported by pre-existing ev¬ 
idence from simulation. Our proposal for a complete 
density profile matching to NFW on large scales is phe¬ 
nomenological, and new to this work. We apply this to 
the dwarf spheroidal (dSph) galaxies, finding limits on 
the axion mass for a solution to the cusp-core prob¬ 
lem. Our choice of data sets and methodology applied 
to the axion core model is entirely new to this work, 
and makes concrete links to cosmological limits. While 
we refer explicitly to our work as on ‘axion DM,’ the re¬ 
sults are more generally applicable to any non-thermal 
scalar DM candidate. We hope that the conclusions we 
draw inspire further study of these models within the 
community. 

We begin in Section with discussion of some in¬ 
tuitive aspects of the soliton profile andthe relation be¬ 
tween the central density, axion mass, and core radius. 
In Section [3.1[ we provide a complete model for the halo 
density profile for axion DM with explicit formulae. We 
use these formulae in Section 13.21 to find limits on the 
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axion mass based on stellar-kinematical data for Fornax 
and Sculptor dSph galaxies taken from WPll. We con¬ 
clude in Section]^ Derivations, pedagogical notes, and 
certain numerical aspects are relegated to the Appendix. 


2 SOLITON CORES FOR SCALAR DM 


Quantum and wave-mechanical properties of axion DM 
allow for pressure support. Bose-Einstein condensation, 
if it occured, could also lead to large correlation lengths. 
Neither phenomenon will occur if the axion is modelled 
as pressureless dust with classical gravitational interac¬ 
tions. Although from the particle point of view, possible 
Bose condensation of the axion field is a quantum phe¬ 
nomenon, from the field point of view it can be studied 
classically | Guth et al.|2014 1. The classical analysis also 
allows one to derive the axion sound speed and resulting 
Jeans scale ( Khlopovet al.| 1985 1. Indeed, the character¬ 
istic wavelength for the Bose-Einstein condensate, the 
scale over which thermalisation occurs and gravitational 
growth is suppressed, is none other than the Jeans scale 
in the fluid description of the classical field under linear 
density perturbations. 


2.1 Schrodinger-Possion System and Ground 
State Solitons 


We work in the non-relativistic approximation, where 
the full Einstein-Klein-Gordon (EKG) equations reduce 
to the Schrodinger-Poisson system ( Seidel fc Suen|19OT 
|Widrow fc Kaiser|1993| ): 


• ft'9V' 


2m.a 


VW = AtxG%1)iP* 


( 1 ) 


The axion mass is mria, and the field i/' is related to the 
WKB amplitude of the axion field, (j) (see Appendix |Al[ ). 

The soliton solutions of this system, and their rel¬ 
ativistic completion, were first studied by |Ruffini "fcl 
[Bonazzola] ( |196^ , who identified the scaling symmetry, 
and the maximum stable mass. High-resolution cosmo¬ 
logical simulations of the Schrodinger picture for DM 
by [Schive et al.| ( |2014a[ ) using an adaptive mesh refine¬ 
ment scheme powered by GPU acceleration have pro¬ 
vided evidence that ultra-light (m^ ~ eV) scalar 

DM halos form kpc scale soliton cores that are stable 
on cosmological time scales. 

Gonsider stationary wave solutions of the form[^ 

( 2 ) 

for a function y(r) which depends on the radial coor¬ 
dinate r, but does not vary with time. Taking drj = 0 


assumes phase coherence and zero fluid velocity. We ex¬ 
pect loss of phase coherence at large distances in real 
systems, and expect that this is related to the transition 
from soliton to NEW profile discussed in SectionjT^ In¬ 
deed, phase coherence is lost outside of the density cores 


in the simulations of Schive et al. (2014a I. 


In spherical co-ordinates we have the following sys¬ 
tem of ODEs: 

2 x' 


X + — = 2(U-7)x, 
r 

V"^ 2 U' 2 

V -f-= X . 


( 3 ) 


where V, x> a-nd 7 are all dimensionless quantities 


defined in Appendix A2 The solutions x(r) are a special 
type of soliton known as an oscillaton (see Appendix [b| 
for clarification of terminology) in the axion field </>, with 
soliton prohle x, and density prohle p = ■ 

In the stationary wave solution, Eq. § , it is im¬ 
portant to note that 7 is a parameter to be solved for, 
related to the total energy. The ground state is the state 
of lowest energy and depends on only one length scale. 
It also has no nodes, and with fixed central density 
this allows us to find the unique value of 7 numerically. 
The stationary wave solution must obey the condition 
Imtpl ^ \htp\ (first order WKB approximation: see Ap¬ 


pendix Al I, which translates to 7 <C 1 in dimensionless 


variables. This must be satished if solutions to Eqs. 
are consistent non-relativistic limits of the solutions to 
the full EKG equations. 

The soliton density profile must have the properties 
p( 0 ) = const, and p{r —>■ 00 ) = 0 , such that p is the solu¬ 
tion of the boundary value problem in Appendix |A5| We 
write the density profile as: 

Psoiir) = f{ar), (4) 

where f{y) has no explicit scales and has the correct 
asymptotic behaviour. Restoring units: 

Psoi{r) = 2m^Mpj/(r/raoi), (5) 

rsoi := -• (6) 

ama 

Solutions to Eqs. possess a scaling symme- 


Jl, 

try ( Ruffini fc Bonazzola] 19691, as discussed in Ap¬ 
pendix M" which is used to fix the appropriate scale 
for astrophysical systems. Upon applying a rescaling by 
A, all dimensionful quantified The density must scale 
A'^p. Rewriting 

Psoi(r 


as p - 


Pcrit 


= 3soif(r/rsoi), 


( 7 ) 


we End the relationship between the soliton density pa¬ 
rameter, (5soi, and its characteristic radius, rsoi: 


<5. = 


5 X 10'‘\ V /i 

0 


h \ / nia Y 
. 7 ) V 10-22 eV/ 


[kpc] 


( 8 ) 


^ Separating the function t) into polar coordinates in 
this way is sometimes referred to as the Madelung transfor¬ 
mation 1 Madelung|1926 l. 


^ Except the axion mass, which does not appear explicitly 
in Eq. §, and hence only sets the units, are affected, so that 
r -I- r/X and rgoi ->■ rgoi/A. 
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The value of a is fixed from numerical fits in Ap¬ 
pendix |A6| 

A rescaling of the density profile is seen to affect 
only (5aoi —> A'^haoi, with the relationship between the 
central density, the scale radius, and the mass com¬ 
pletely being fixed. The numerical value of a depends 
on the choice of the functional form of f{ar). Choos¬ 
ing a different definition for raoi, e.g. the half-density 
radius used by Schive et al. ( 2014a|b' l, or even a com¬ 
pletely different functional form for / will only change 
the numerical co-efficient in Eq. ([^ and not the func¬ 
tional relationship between 5aoi and raoi, or their depen¬ 
dence on the axion mass. These features are universal. 


In Appendix A5 we fit the form of f{ar) from numerical 


solutions, and fix the value of a for our chosen fit. 

The axion mass can be re-expressed in terms of the 
linear Jeans scale, ryun/kpc oc (mo/10“^^ (Hu 

et al. 2000j ). Substituting into Eqs. ( |7|8| ), the scale radius 
of any soliton with central density psoi(O) is given by 


r-aol OC - ryiin . 

V Pcrit / 


( 9 ) 


This simple result fixes the scaling properties of axion 
density cores using only the scaling symmetry of the 
SP system. The existence of the scale in the ground 
state was guaranteed by the uncertainty principle, which 
holds by virtue of the large occupation numbers and 
the classical wave-mechanical description of the axion 
field. The scaling symmetry occurs because in the non- 
relativistic limit there is no scale in the SP system. 
Such a scaling should therefore hold for all large occupa¬ 
tion number, non-relativistic axion/scalar DM density 
configurations in the small-radius and long-time limits. 
That is, soliton cores with this scaling property are uni¬ 
versal features of such DM models. 


2.2 Soliton Core Radius 

We now consider the relation of the soliton radius to the 
de Broglie wavelength. Consider a particle of mass nia 
on a circular orbit around the soliton core. Its velocity 
is given by: 


2 GM{< r) 

r 

(10) 

The uncertainty principle requires: 


pr ^ h 

(11) 

The de Broglie scale is the point where prun 
thus: 

= h, and 

IGM{< ruB) , 

rdBrUaX — = h 

V rdB 

(12) 

1 

^ M{< rds) = 2 ^ 

miG rdB 

(13) 


In Fig.[^we solve this equation graphically for some typ¬ 
ical dwarf galaxy parameters and find ran = 1.04 kpc. 
The de Broglie wavelength found in this way is close 



Figure 1. Solving Eq. | |13[ l. Blue: left hand side, the mass 
of the soliton as a function of radius for some representative 
dwarf galaxy parameters. Green: right hand side. The inter¬ 
section solves for the de Broglie radius thus defined. The red 
dot shows the core radius as defined by |Schive et al.| | |2014a^ , 
ri/ 2 , while the green dot shows the de Broglie radius ob¬ 
tained from the uncertainty principle. Due to the steeply 
falling soliton profile outside of the core, neither definition 
of core radius fixes the density of the point of transition to 
NEW accurately l|Schive et al.|2014a|b[|. 


to the core radius at half central density, ri/ 2 , where 
Psoi(ri/2) — Psoi(0)/2. 

The form of the scaling in Eq. is the same as 
for the halo Jeans scale extrapolated from linear theory 


( Hu et a'L]|2000 Marsh fc Silk|[20l3 Guth et al.||2014 i. 
Both follow from the scale symmetry and dimensional 
analysis, and thus their interpretation as related to the 


de Broglie scale remains the same (see also Hlozek et al. 


20141. Beyond scaling alone, how exactly are the lin¬ 


ear Jeans scale and the soliton size related numerically? 
Can this be used to predict properties of axion halos, 
including the core-size/mass relation, and the transition 
to NEW? 


The simulation results of Schive et al. (2014a I in¬ 


dicate a soliton core connecting to an NEW profile at 
larger radii. No prescription for the matching radius is 


provided, however. In a follow up paper (Schive et al. 
|2014b[ ), the same authors assessed halo formation in 
this model as the gravitational collapse of a system of 
solitons. Using the uncertainty principle and the virial 
theorem they motivated a relationship between the to¬ 
tal halo mass, Mh, the soliton core mass M{< ^ 1 / 2 ), 
and the minimum halo mass, Aimin'- AI{< ri/ 2 ) oc 
{Mh/MmiA^Mm'm- _ 


Hu et al. (20001 and Marsh & Silk (20131 related 


the core radius to the halo Jeans scale. Can the same 
quantity be used to define a matching point between 
soliton and NEW? With respect to the linear Jeans scale 
rj, the halo Jeans scale, rj^h, satisfies 


P(rj,h, 

Pcrit 


rj 

rj,h 


(14) 
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dSph 

logio(o-^/km^ s 

Error 

logic (r/kpc) 

Error 

Fornax 

2.00 

0.05 

-0.26 

0.04 


2.32 

0.04 

-0.05 

0.04 

Sculptor 

1.62 

0.06 

-0.78 

0.04 


2.13 

0.05 

-0.52 

0.04 


Table 1. Data used in this work, taken from WPll. We approximate the likelihoods to be two-dimensional uncorrelated 
Gaussians in log]^Q(r, cr) for each data point. 


In order for this relation to define a matching point 
between soliton and NFW profiles, both profiles must 
provide a solution to this equation at the same radius. 
The logarithmic slope, T, of the NFW profile is in the 
range — 3 < F < — 1 and so there is always a single 
unique solution for rj^h- However, the slope of the soli¬ 
ton profile is zero at the origin, decreasing to F < —4 at 
large radius. For the soliton, there can be either zero, 
one or two solutions for ryj,, with a single unique solu¬ 
tion at the point where F = —4. Using the relationship 
in Eq. ^ it is possible to show that solutions only occur 
for Ss ^ 1. Halos, on the other hand, correspond to non¬ 
linear density perturbations with 5s ^ 1. Therefore, the 
transition from soliton to NFW profile in a halo must 
occur on radii smaller than the linear Jeans scale extrap¬ 
olated using the local density, i.e. on r < rj,h defined by 
Eq. (141. This implies that soliton cores are more com¬ 
pact than expected from linear theory. This will turn out 
to have important implications for a possible solution 
to the cusp-core problem using ULAs/scalar field DM. 


3 AXION HALO DENSITY PROFILES AND 
DWARF GALAXY CORES 

In this section we will first define our proposal for the 
complete halo density profile of axion /scalar-field DM. 
We then go on to estimate the parameters in this pro¬ 
file using the Fornax and Sculptor dSph density profile 
slopes as measured by WPll. We introduce this mea¬ 
surement and the approximations it uses, define a like¬ 
lihood from this, and perform a Monte Carlo Markov 
chain (MCMC) analysis. The results are all contained 
in Table and Figs. and 


3.1 Defining the density profile 


The halo density profiles for ultra light scalar dark mat¬ 
ter observed in the simulations of |Schive et al.| | |2014a|b] ) 
consist of an NFW-like outer region, with a prominent 
solitonic core. The observed transition is sharp, and we 
model it as a step function: 

p(r) = e(re - r)psoi(r) -I- 0(r - re)pNFw(r) . (15) 


The transition to an NFW profile at some radius is 
expected on a number of physical grounds. Axion/scalar 
field DM is indistinguishable from CDM on large scales, 
and this is the basis of the Schrodinger approach to sim¬ 
ulating CDM (Widrow fc Kaiser|1993 Uhlemann et al. 


20141. Thus, on scales much larger than the de Broglie 


wavelength, the scalar field halos should resemble those 
found in N-body simulations, i.e. NFW. Furthermore, 
no long range correlation should occur; the smallest ob¬ 
jects are solitons, with a characteristic granularity fixed 
by the Jeans scale (Schive et al. 2014a|b Guth et al. 


20141. On large scales, phase decoherence should oc¬ 


cur, violating the assumption in our soliton ansatz that 

dr'y = 0 . 

The dynamics will be that of an interacting gas 
of solitons in a decoherent scalar field background. On 
length scales larger than the soliton radius the usual 


thermodynamic arguments relating to dust (Binney & 


Tremaine 20081 will apply, leading to an NFW-like 


profile (on large scales this equivalence between the 
Schrodinger picture and dust thermodynamics can be 
derived using the Wigner distribution,e.g. [Widrow fc| 
Kaiser 1993 Uhlemann et al.]|2014 (. The transition in 
behaviour from soliton to NFW governed by the deco¬ 
herence scale, and the mass function of solitons in the 
outer halo are all interesting questions, and will be the 
subjects of forthcoming papers. The NFW density pro¬ 
file is given by ( Navarro et al.|1997 l 

PNFw(r) _ S 

= {r/rAA + r/rA^ ' ^ 

We fit for the soliton density profile in Appendix |A6| 

Paol(^) _ __ ('ir.N 

Pcrit (1 + (r/rsoi)h« ■ 

We take Sa as a free parameter and rgoi is fixed in terms 
of it by Eq. ^ using a = 0.230. The rescaling param¬ 
eter, A, can then be found by setting Csoi = {XamA~^■ 
We recall that consistent solutions with | 7 /m| <C 1 re¬ 
quire A < 1, which also guarantees <j>{0) < O.SMpi and 
the stability of the soliton ( Seidel fc Suen[]1990 (. 

We match the soliton and NFW profiles at a fixed 
value of the overdensity, S = eSa, defining the matching 
radius r^. Since we have not been able to find an accu¬ 
rate estimate that fixes e analytically, we take it as a free 
parameter [^Continuity of the density at this point fixes 
one of the two free parameters in the NFW profile. We 
choose this to be Jchar and take the scale radius to be a 
free parameter. For a given ULA mass, our halo density 
profile thus has three additional free parameters: 


{<53,e,)-4 . 


(18) 


^ An order of magnitude estimate for the matching radius 
based on the de Broglie scale will not be enough. The soliton 
density falls rapidly for r > ^rid so 0(1) numerical co¬ 
efficients have a large effect on the estimate for e. 
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Parameter 

Prior 

logio(m<i/eV) 

U{X,-19) 

loglO 

f/(0,10) 

logio e'’’® 

t/(-5,logio0.5) 

loglo('’s’®/kpc) 

t/(-l,2) 


Table 2. Priors on our density profile parameters, where 
U{a, b) is the uniform distribution on [a, b]. The lower bound 
on the mass prior, X is given two different values: Xcmb = 
—25 plozek et al.|p014^ and XjjuDF = ~23 l |Bozek et al.| 
|2014|l, referred to as CMB and HUDF priors respectively. 


Parameter 

Posterior 

(ma/eVjnuDF 
(ma/eV)cMB 
logic <5? 
logic 

< 1.1 X 10“^^ 

< 1.0 X 10-22 

(95% C.L.) 
(95% C.L.) 

c cc:+0.06 

0.00_Q Qg 

6.19 ±0.05 


Table 3. Posteriors on the constrained density profile pa¬ 
rameters. The mass constraint is quoted for both the CMB 
and HUDF priors (see Table |^. Upper and lower errors on 
the central densities are given as the 16th and 84th per¬ 
centiles and are quoted for the HUDF prior only (see text 
for discussion). The upper bound on the axion mass reflects 
the minimum core size consistent with observations. 


We emphasise again that a theoretical model for 
the value of e, which could depend on redshift, central 
density and/or particle mass, could perhaps be derived. 
In this case, the density profile has just as many free 
parameters as an NFW prohle. Thus, in any given halo, 
a core measurement would predict the point of tran¬ 
sition to NFW, and a measurement of the outer halo 
fixing concentration and scale radius would predict the 
corresponding inner core size. 

The data we use imply cored profiles on the ob¬ 
served radii, and so in our phenomenological model with 
free e we will find that the NFW parameters e and are 
unconstrained. This shows that the data prefer soliton 
cores to NFW prohles. For a complete Bayesian analysis, 
we include and marginalise over the NFW parameters 
in our constraints, as described in the next subsection. 


3.2 Fitting to Fornax and Sculptor 


While there are several ways of analysing the observed 
data, we will focus on the method used by WPll who 
measured the slopes of dSph mass profiles directly from 
stellar spectroscopic data. They use the fact that some 
dSphs have been shown to have at least two stellar pop¬ 
ulations that are chemo-dynamically distinct (see e.g. 


Tolstoy et al.|2004{|Battaglia et al.|2006|[Battaglia et al.| 


20111. Measuring the halflight radii and velocity disper¬ 


sions of two such populations in a dSph allows one to 
infer the slope of the mass profile. The method of WPll 
has the advantage that it does not need to adopt any 
dark matter halo model a priori. Therefore the data can 
be used to test theoretical density profile models with¬ 
out having to run computationally expensive hts to the 


full stellar data for each value of the theoretical param¬ 
eters. 

We model the results of WPll as providing two- 
dimensional Gaussian distributions for the half light 
radii, rh,i, and velocity dispersions, for each of 

the two stellar populations, i, in Fornax and Sculptor. 
While the exact results display some covariance between 
o and Th the Gaussian approximation is far simpler to 
analyse, and accurate enough for the purposes of this 
study. This follows the approach taken by |Lombriser fc| 
Pefiarrubia (20141 for testing chameleon gravity using 
dSphs. The data we use are given in Table 

For our density profile the mass internal to any ra¬ 
dius, M(< r), can be computed analytically (and so can 
the derivative of the mass, dM/dr). These are the only 
ingredients necessary in anlaysing this simplified version 
of the stellar kinematic data. The velocity dispersion at 
the stellar half-light radius obeys the following empirical 
relationship: 


= 


2GM{< rn) 
5rh 


(19) 


This relationship is related to the virial theorem, and is 
found by solving the Jeans equation and finding a ‘sweet 
spot’ where a wide variety of density and velocity pro¬ 
files agree (including projection and anisotropy effects). 
In this work we do not perform a full Jeans analysis 
and use only this analytic relationship and the derived 
errors on it, following WPll. 

The two-dimensional nature of the (r, a) data can 
be accounted for using an effective one-dimensional er¬ 
ror (e.g. Ma et ar]|2013 1. If the data has central values 
(x, y) and standard deviations (Sx, ^y) then for a given 
model y = f{x) the effective one-dimensional error is 
given by: 


Sefl — + 


dfjx) 

dx 


The likelihood, £., can then be approximated as 

'-{fix)-yf' 


C oc exp 




( 20 ) 


( 21 ) 


We use the data for Fornax and Sculptor with equal 
weight in the combined likelihood. The axion mass, ma, 
is a global parameter which is the same for both Fornax 
and Sculptor. Our final model thus has seven parameters 
in total: 


r fF jtS F S S Ft 
\ma, Og, 0 g,Tg,Tg,€ , e J , 


( 22 ) 


with F, S, labelling Fornax and Sculptor respectively. 

The priors on these parameters are given in Table[^ 
We take Jeffreys’, or ‘least information’, priors in a fixed 
range for all parameters. We will discuss our results in 
detail shortly, but mention here those aspects relevant to 
priors. We find that the NFW parameters (e, r^) are un¬ 
constrained. The priors on these parameters are there¬ 
fore irrelevant in quoting marginalised constraints on 
other parameters and we take them extremely wide to 
explore all possibilities. We impose an upper bound on 
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Figure 2. Two and one-dimensional marginalised posteriors on the constrained density profile parameters for the HUDF 
mass prior (see Table [^. Contour levels show [0.5,1.0,1.5, 2.Ojc preferred regions. Vertical dashed lines show [0.05,0.5,0.95] 
percentiles. This plot is made using triangleplot (|Foreman-Mackey|2014[|. 


the e prior of e = 0.5, so that the match must occur 
outside the half-density radius, consistent with the sim- 


bound and relies only on linear constraints from the 
CMB. It is thus extremely reliable. Our alternative prior 


ulation results of Schive et al. (2014a I. As already men- uses the results Bozek et al. (20141, which constrains 


tioned, the fact that these parameters are unconstrained 
has physical meaning: our MCMC finds no peak corre¬ 
sponding to an NFW profile, and thus the data prefer 
soliton cores to NFW cusps. The central densities are 
well constrained and so the results are independent of 
the prior. 

We find a one-sided constraint on the axion mass, 
and so the percentiles quoted depend on the prior for 
the lower bound. We take two such priors. The first uses 


the results of Hlozek et al. (2014l, which place an ap¬ 


proximate lower bound on ULAs to be all of the DM 
of TUa > 10“^® eV. This is a very conservative lower 


rria > 10“^® eV at more than 8 <t significance using 
Hubble Ultra-Deep-Field (HUDF). While this is a very 
strong bound, it relies on more assumptions about the 
astrophysics of reionisation and on the non-linear struc¬ 
ture formation of ULAs. This HUDF bound is, however, 
still rather conservative compared to other constraints 
to ULAs using non-linear scales (e.g. Lyman-alpha for- 
Amendola fc Barbieri|2006 1 . 


est 


We evaluate the likelihood using emcee, an afBne- 


invariant MCMC ensemble sampler (Foreman-Mackey 
et al.|20l3 1. We use 200 ‘walkers.’ Convergence is tested 
by evaluating the auto-correlation time, which is found 
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Figure 3. Velocity dispersions are estimated at the stellar half-light radii from M(< r^), Eq. | |19[ l. We show velocity dispersions 
and density profiles of 200 random samples from our MCMC using the HUDF mass prior (see Table[^. The large core in Fornax 
favours low axion masses. Consistency with HUDF predicts a rapid turnover in the density slope at radii slightly larger than 
those observed, caused by the transition from soliton to NFW. Data from WPll, quoted in Table[^ 


to be 0(80) for each parameter. With 5000 steps per 
walker the likelihood is thus evaluated over many cor¬ 
relation times, giving many independent samples. We 


have used an ‘MCMC hammer’ (Foreman-Mackey et al. 


20131 to crack a very simple parameter-constraint-nut, 


but this gives us a high degree of confidence in our re¬ 
sults, and allows us to present them in a completely 
Bayesian fashion. 


Our results are presented in Table and shown for 
the HUDF priors in Fig. As already mentioned, the 
NFW parameters {rs,e) are unconstrained, and so we 
do not show them. A weak constraint could be found on 
these parameters if one were to solve the full Jeans equa¬ 
tion to fit the stellar kinematic data, rather than using 
only the empirical relationship Eq. (191. The constraint 
is caused by requiring a fixed enclosed mass at large 
radius. This consequently also imposes a weak lower 
bound on the axion mass to give a finite core radius, 
as in Schive et al. (2014aI. Such an analysis is beyond 


the scope of this paper. Information about the total halo 
virial mass could also allow one to impose the core-halo 


relation of Schive et al. (2014b I, introducing a degener¬ 
acy between and e. 


For the central densities, shifts in the central val¬ 
ues using the CMB prior compared to the HUDF prior 
are smaller than la for Fornax and Sculptor, with com¬ 
parable errors in both cases. The shift for Fornax is 
ever-so-slightly larger due to the degeneracy between 
6s and rria enforced by Eq. ([^ having an effect when 
marginalising over rUa with a different prior in the un¬ 
constrained, low mass region. The degeneracy between 
these parameters is more pronounced for Fornax as it is 
less cored than Sculptor. This same effect is the cause 
of the non-Gaussianity in the one-dimensional central 
density distribution for Fornax. 


Random samples from our emcee runs showing the 
density and velocity profiles colour-coded by axion mass 
are shown in Fig. for the HUDF prior. This makes 
one clear point: since the HUDF prior already restricts 
the mass, cores cannot be too large. Therefore, if ULAs 
are responsible for dSph cores, velocity dispersions and 
density profiles must turn over quite rapidly on radii 
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just larger than those observed. This is, in principle, a 
falsihable prediction of the ULA model for cores. 


4 DISCUSSION AND CONCLUSIONS 


In this paper we have investigated stable solitonic den¬ 
sity profiles, which are expected to be the smallest struc¬ 
tures formed in models of axion/scalar field DM. We 
found these objects via numerical solution of the bound¬ 
ary value problem in the Schrodinger approach to DM. 
Via use of a scaling symmetry we found a universal re¬ 
lationship between the central density of the soliton, 
the soliton radius and the DM particle mass. We fur¬ 
ther explained the soliton size by recourse to the un¬ 
certainty principle. The standard thermodynamic argu¬ 
ments for dust, and the formal equivalence between the 
Schrodinger picture and dust on large scales, further 
supported by evidence from simulation, lead us to con¬ 
sider solitons as the central cores in NFW halos. We pre¬ 
sented a new phenomenological model for such a density 
profile. 

If the DM is ultra-light, then solitons may be re¬ 
sponsible for the density cores observed in dSph galax¬ 
ies, with soliton density profiles in the inner regions, 
and an outer NFW profile. We investigated the valid¬ 
ity of this claim by performing an MCMC analysis of 


the simplified stellar kinematic data of WPll (Walker 
fc Penarrubia|20TT I for Fornax and Sculptor. This data 


shows a preference for cores over cusps, but other stud¬ 


ies prefer cusps (e.g. Breddels & Helmi 2013 Richardson 
fc Fairbairn|2014 l. The preference of the WPll data for 


large cores shows up as an upper bound on the axion 
DM mass of nia < 1.1 x 10“^^ eV at 95% C.L., leaving 
the NFW parameters unconstrained and marginalised 
over. The axion mass bound applies only if the dSph 
cores are solitonic, and not if they are caused by bary- 
onic effects for standard CDM (be it composed of heav¬ 
ier axions, WIMPs etc.). This bound is consistent with 


the best ht mass m, 


= 8.111;? X 10" 


eV (an approx¬ 
imate 95% C.L. region 0.47 < ma/10“^^ eV < 1.13) 


of Schive et al. (2014a I found by solving a simplified 


version of the full Jeans equations for a single stellar 


sub-component in Fornax alone. Schive et al. 12014a) 


also checked the (rough) consistency of their results with 
density profile slope measurements, but our analysis is 
the first to use these data together in a complete and 
Bayesian manner. In future we hope to conduct a full 
Jeans analysis for our model, following the similar anal¬ 
ysis for self-interacting scalar field DM by |Diez-Tejedor| 

[eraI1 ( [200T4| ). 

The existence of an upper bound on the DM par¬ 
ticle mass in the axion model has a number of interest¬ 
ing consequences. Structure formation in this model is 
substantially different from CDM due to the cut-off in 
the power spectrum on small scales caused by the ax¬ 
ion sound speed and resulting Jeans scale. |Bozek et al.| 
(2014) found that, if the DM is composed entirely of 


ultra-light axions, rua = 10 eV is just consistent 
with HUDF. Our upper bound is right on the edge of 


this limit, and suggests two possible outcomes warrant¬ 
ing further study: 


• Pessimistic: Structure formation and the require¬ 
ment of dSph cores put conflicting demands on axion 
DM. This places the model in a ‘Catch 22’ analogous 
to WDM^ This particle physics model is no longer a 
catch-all solution to the small-scale crises, and addi¬ 
tional mechanisms are required for its consistency. 

• Optimistic: Ultra-light axions are responsible for 
dSph cores. The cut-off in the power spectrum is just 
outside current observational reach. Near-future experi¬ 
ments will turn up striking evidence for axions in struc¬ 
ture formation and in study of the high-z universe. 


In the pessimistic case, axions have nothing to do 
with core formation in dSph galaxies, and another mech¬ 
anism is needed to form the cores. The DM can still 
be composed of axions, but there is now no limit on 
the mass based on the dSph observations. Cores may 
be formed by baryonic processes, putting axion DM 
in the same situation as any other DM model, e.g. 
WIMPs. Cores could also be formed by adding strong 
self-interactions for CDM or for axions. One could also 
modify gravity, for example if the modulus partner of 


the axion played the role of a chameleon field (Khoury & 


Weltman|2004 Lombriser fc Penarrubia|2014 1. Another 

alternative would be to keep cores from axion solitons, 
but compensate structure formation by a boost in the 
primordial power or some other alteration to cosmology. 
A compensation based on a mixed DM model, however, 
is unlikely to succeed ( [Marsh &: Silk|2013p . 

There are still opportunities to discover evidence 
for axion DM in the pessimistic case. We briefly out¬ 
line some of these for the ‘most pessimistic’ case (from 
a particle physics view point) that cores are formed en¬ 
tirely by baryonic processes, and consider only searches 
influenced in some way by the Jeans scale or soli¬ 
tons, i.e. by non-trivial gravitational behaviour. For 
10“^^ eV < rria < 10“^° eV axions are consistent with 
current constraints from cosmology, but the cut off in 
the power spectrum may still be observable using mea¬ 
surements of the weak lensing or 21cm power spec¬ 
tra. Axions can further be evidenced by their effect 
on black holes via the superradiant instability ( |Arvan-| 
itaki fc Dubovsk^|2011 1 . Spinning super-massive black 
holes indirectly probe/constrain masses in the range 


10 " 


eV < rria 10 “ eV (Pani et al. 


2012 


Brito 


et al. 20141. More direct signatures in gravitational wave 


detection may come from solar mass black hole s for ax- 
ions in the range 10“^® eV < ma < 10“^° eV (Arvani- 
taki et ^[2014 ). 


Axion DM of any mass will contain solitons some¬ 
where in the mass spectrum at some point in cosmic 
history. This could have a number of consequences rel¬ 
evant to both the optimistic and pessimistic scenarios 


Phrase coined by [Maccib et al.j ( [2012^ . For ULAs, one 
might call it a Catch 10“^^. 
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outlined above. This model predicts small and dense 
cores in all DM halos, with core size inversely related 
to central density and halo mass (Eq. and Schive 


et al. 2014bI. Snch cores may provide seeds for high- 


2 quasars (Schive et al. 2014bI. The granular nature 


of DM composed of solitons could be detected by mea¬ 
sures of substructure (which we discuss further below). 
DM composed entirely of solitons/oscillons/boson stars 
may have a number of interesting features for gravita¬ 
tional wave observations, as discussed in, e.g. [Evslin ^ 
|Gudnas^ ( |2012] ); |Macedo et al.| ( |2013[ ), that distinguish 
soliton DM from an equivalent model with black holes, 
due to the absence of a horizon. The soliton content of 
an axion DM halo will also impact the direct detection 


prospects (Hoskins et al. 20111, while the evaporation 


of solitons from axion self interactions enhanced by the 
high number densities could have indirect signatures. 
Other signatures of the soliton components of an ax¬ 
ion halo, e.g. in precision time-delay experiments with 
atomic clocks, may be similar to those with topological 
defect DM (Pospelov et al|2013 Derevianko & Pospelov 


2014 Stadnik fc Flambaum|20T4 1. 


There are also ‘intermediate’ cases to consider. 
ULAs with rua < 10“^® eV could be detected as a 
sub-dominant component of the DM with high-precision 
studies ( [Hlozek et al.|[2014[ ). A sub-dominant compo¬ 
nent at these masses is likely the expectation based 
on string models with sub-Planckian axion decay con¬ 
stants (although see [Bachlechner et al.||201^ . ULAs 
with TUa > 10~^^ eV as the dominant DM could also 
play some role in core formation alongside inefficient 
baryonic processes. 

Finally, we turn to discussion of the optimistic case, 
where axions at the edge of our allowed region in mass 
are solely responsible for dSph cores. In this case, soli¬ 
ton cores have a large mass 0(10® Mq) and large radial 
extent 0(1 kpc), and the cut off in the halo mass func¬ 
tion (in the field) occurs at M ~ 10® Mq (the 2 = 0 
Jeans mass, Marsh fc Silk|20T3 Bozek et al.|2014 |. 

We have demonstrated the ability of soliton cores 
to fit the dSph observations, and suggested a complete 
density profile including the NEW piece. This density 
profile could be fit to many other observations, includ¬ 
ing rotation curves (de Blok et al.|20'oT Robles & Matos 


20121 and stellar clump survival (Lora fc Magana|20 141. 


We argue that further study of this model should use 
the density profiles we advocate in this work, and not 
‘multi-state’ or other adaptations of the scalar field 
DM model. The outer NEW piece obviates the need 
for such additions, and is consistent with many lines 
of reasoning. Further study should, as we have in this 
work, account for consistency of the cosmology, as well 
as between different density profile constraints. Consis¬ 
tency with cosmology and our core size fits in this work 
suggest that density profiles just outside 1 kpc should 
be much steeper than cores and make a transition to 
NEW shortly thereafter. It should also be investigated 
whether the scale symmetry for soliton cores can be used 


to explain the observed scaling relationships in dwarf 
galaxies ( Kormendy fc Freeman|[2014 Burkert|2014 1. 

While density profiles provide the motivation for 
rUa ~ 10“^^ eV, cosmology suggests that the most strik¬ 
ing evidence for such a model will be found in sub¬ 
structure and high -2 galaxy formation. Measures of sub¬ 
structure from mili-lensing ([Dalai fc Kochanek [2002 I, 
tidal tails (Johnston et al. [2002 1, and strong lensing 
( Hezaveh et al.|2014 l should be sensitive to the absence 
of small-scale structure predicted by an axion solution to 
the cusp-core problem. Further study is required to for¬ 
mulate the predictions of the axion model in this regard. 
This should come from simulation, but also from semi- 
analytic study of the conditional mass function and ex¬ 
tended Press-Schecter formalism ( Lacey fc Cole[[l993 l. 
This study will further address the questions of whether 
axion DM can provide a consistent solution to other 
small-scale crises like ‘missing satellites.’ 

A final set of ‘smoking gun’ signatures in the op- 
timist i c scen ario are suggested by the results of [Bozel^ 


et al. 


(20141. Since nia ~ 10 eV is just consistent 


with the HUDF UV luminosity function, and also just 
able to provide a dSph core, this makes the prediction 
that a JWST measurement of the UV luminosity func¬ 
tion sensitive to fainter magnitudes will only see a value 
consistent with HUDF, and not the larger value pre¬ 
dicted by CDM. The cut off in the mass function also 
makes predictions about reionisation that are less sen¬ 
sitive to astrophysical uncertainties than in CDM. An 
axion solution to the cusp-core problem predicts that 


CMB polarisation experiments ( Calabrese et al. 20141 
should measure: a low value of r ~ O.OSj^a low redshift 
of reionisation 2re ~ 7, and thin width of reionisation 
Szye ~ 1.5. These effects on reionisation, caused by the 
rapid build up of massive structures at the mass func¬ 
tion cut-off, should also turn up striking effects in the 
21cm power spectrum (building on Kadota et al.[[2014 


[Sitwell et al.[[20T^ . Further study of all of these signa¬ 
tures is underway. 
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APPENDIX A: THE 
SCHRODINGER-POISSON SYSTEM 

A1 Derivation 

We begin by deriving our workhorse system, the 
Schodinger-Poission (SP) form of the field equations, 


following Widrow & Kaiser (19931. Further reading on 
this formalism can be found in [Coles (|2003|); |Coles fc] 


icer| ( 

ihig 


Spenced (20031; Johnston et al. (20091, and references 


thereinljThe SP system can be viewed as a fundamental 
picture of structure formation for axions or other scalar 
DM, or as a numerical procedure to model and study 
CDM with a suitable cut-off, having certain advantages 
over N-body simulations ( |Uhlemann et ah ]M4| ). 

The Klein-Gordon (KG) equation for a free homo¬ 
geneous scalar held ij) of mass nia, in an expanding uni¬ 
verse with Hubble rate H, is: 


0 -f 3R<(. -b 


h? 


1 = 0. 


(Al) 


From this point onwards, we will set c = 1, but keep 
factors of h for now. In the regime H m (valid at 
late times, when the axion held is oscillating and be¬ 
having as DM) the equation above can be solved by 
WKB methods: 


h 


V^ma 


(V’e 




(A2) 


for the function ip that is slowly varying with time, i.e. 
\mip\ » \hip\. 

The time-time metric component is (in the Newto¬ 
nian gauge and in physical time): 


5oo = -[l + 2K(r)], 


(A3) 


where V{r) is the Newtonian potential. Treating this 
perturbatively, the KG equation for the inhomogeneous 
held becomes: 






h? 


— 0, 


0 -b ?,Htp - d^di<p -b (1 -b 2V)^(P = 0 . 


(A4) 


(AS) 


In the non-relativistic limit, this equation has the same 
ansatz solution as the homogeneous held equation. Plug¬ 
ging in our ansatz and neglecting terms of order 0{ip), 


® Here, we are using a Schrodinger system as a change of 
variables to understand a classical wave equation, which also 
has a fluid description. Similarly, an understanding of clas¬ 
sical fluids can help in understanding aspects of quantum 


mechanics. See Park (1979 i for further discussion. 


since ip is a, slowly varying function of time, we hnd 
Ip =~j2i -b 




V^<6= ^ 


y2r: 


+. (A7) 


In the regime H m, the second term in Eq. (A51 can 
be neglected. This gives: 


\/2i {—Ip e 


‘l —imatjh 


+ip*e‘ 


irriat / h\ 


V2n 


(vV 






e ““'“1=0. (A8) 


Multiplying the equation above by hj^/2 and identifying 
all terms with the same exponential factor, we hnd the 
familiar Schrodinger equation: 


ihip = --—V Ip + rUaVip . 
2m 


(A9) 


In the sub-horizon, non-relativistic limit, in any gauge, 
the Poisson equation is: 


vV = AtvGp = AnGipip* . 


(AlO) 


Eqs. (A9l and (AlOI form the coupled SP system for 
inhomogeneous scalar helds, give as Eq. ([^. The non- 
relativistic limit in this case caused us to drop the gra¬ 
dient energy from the right hand side of the Poisson 
equation, setting p = \ip\^■ This is valid in the limit that 
k/ma ^ 1 for wavenumber fc, i.e. dx(p/(Mpima) <C 1. 
Consistent solutions to the SP system must respect this 
limit. 

In deriving this form of the held equations, the only 
quantity that has been treated perturbatively is the po¬ 
tential, V, i.e. we are working in the weak-held, Newto¬ 
nian limit of General Relativity. We have not treated the 
held, (p, or the density, p, perturbatively, and therefore 
our results will be valid in the non-linear (in terms of 
density fluctuations about the critical density) regime 
of gravitational collapse. The only clear limitation of 
our solutions will be instability to black hole formation: 
the Jeans limit (for boson stars, this is the analogue of 
the Chandrasekar limit, see [RufHrh_j£^oirnzzolaJ 1969 


[Khlopov et al. 1 1985] I Seidel fc Suen|1990||Choptuik 1993 

for more details). 

The SP system is simply another way of rewrit¬ 
ing the KG equation, and is related to the more famil¬ 
iar fluid treatment of scalar helds in cosmology (e.g. 
Hu||1998 (. The standard fluid picture is most useful for 


studying linear perturbations in Fourier space via the 


sound speed and Jeans instability (e.g. Hwang & Noh 
2009| |Noh et alT]|2013| |Hlozek et al.||2(114| [Alcubierre 


et al.|2015| ).' We will hnd the Schrodinger picture to be 


more intuitive and useful for applications in real space 
relevant to non-linear densities in halos. 
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A2 Dimensionless Variables 


Plugging the soliton ansatz, Eq. it. into the SP system, 
Eqs. 0 . yields: 


v72 2 m<, 7 

)x, 

vV = AtvGx^ ■ 


(All) 


From Eq. ( |A2[ ), we know the dimension of [x] = ['*/’] = 
[ma][<))] = [rria]^, and we will now replace Xi 7 by 

the dimensionless variables Xi ^ s^nd 7 : 

r fh/niaC, 

7—>-7ma/fi. (A12) 


In these dimensionless variables the SP system is: 

V'x = 2(F-7)x, (A13) 

VV = x" . (A14) 


Tildes can then be dropped. 


A3 Scaling Relations 

The SP system of Eqs. obeys a scaling relation 
(r, X, P, 7 ) —> {r/X,X^xA for scale factor A 

( RufBni fc Bonazzola||196^ . This is easy to check: 

w2 1 vv2 1 ~ 

V X = X = X, 

2(V-7)x = 2 ^^^J, 

V\ = 2{V-x)x, 

x= ^(fo- 7 )x, 

77 V^x = 2(V-7)x. 

We can also extend this scaling relation to include the 
mass of the soliton: 

rrs pTs 

Ms = I AnPpdr = I dvrr^lxl^dr 

= (A15) 

and the density of the soliton: 

Psoi = Ixl^ = ^Ixl^ = ^ (A16) 

Therefore, we can summarize the scaling relation as: 

(r, X, R, 7 , Ms,ps) ^ (r/A, A"x, A V, A" 7 , AM,, A^p.) 

(A17) 

This scaling symmetry is very powerful, and makes 
results about the small-scale clustering of axion/scalar- 
field DM have a universal character. The scaling can 
be used to find the density profiles on different length 
and density scales at fixed axion mass, rria, but can also 
be used to rescale the mass at fixed length scale. The 
scaling symmetry will be useful to us in Appendix | A5| 


where we numerically find soliton solutions with arbi¬ 
trary central density p and then scale them to astro- 
physical densities. 

For the SP system to be used as a model for the 
gravitational collapse of non-relativistic DM axions, its 
solutions must be consistent limits of the full KG equa¬ 
tion satisfying 7^1 and k/m <C 1. Solutions to the 
SP system for a Hxed scaling are formally solutions re¬ 
gardless of the value of 7 and need not satisfy these 
constraints. However, when we fix the scaling for real 
astrophysical systems, the rescaled values must satisfy 
these constraint. In typical DM halos, the virial velocity 
is u ~ 100 km s“^ ^ c. Since both 7 and k are related 
to the kinetic energy, we therefore expect to find these 
limits satisfied in astrophysical objects. In our explicit 
examples of density profiles for dSphs this is indeed the 
case. 


A4 Power Law Solutions 


In this subsection, we are especially interested in the os- 
cillaton profile at small radii. In the limit r —>■ 0, we can 
assume that the dominant term in the series expansion 
of the solution xC^") is of the form x = Cr ^, for a given 
constant C and exponent j3 to be determined below. 

In the case jS yf {—1,0} the Schrodinger equation 
implies: 

x" + 2 ^ =/3(^-kl)C’r'’“", 
r 

= 2{V - x)Cr^ ■ (A18) 

Therefore, as r —>■ 0, we find that V (r) obeys the equa¬ 
tion: 


^ 2 r 2 


(A19) 


Substituting into the Poisson equation: 


+ i-) 

r4 


V +2— = Cr 
r 

,(/?+!) 

r-4 


- 2/3 ^^ { 


PiP +1) = CA 


- /^2^2/3-|-4 


we find that the SP system has a solution for P = —2. 
This implies a steeply increasing central density prohle: 

P{r) = R) ~ (A20) 


However, in this case, we can see that the mass within 
any radius r, diverges: 


M = 


rr^ 

/ p(r)47rr" 

Jo 


dr ~ 47r I — 


(A21) 


Therefore, the requirement that the central mass is finite 
excludes the case P = —2 . 

Next we check the cases /3 = 0 and P = —1 sepa¬ 
rately. 

For P = —1, we have that x = C/r, x! = ~C/R, 
xP' — 2C/r®, and thus: 




0 = 2(V-7)-. 

r 
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This expression must be true as r —>■ 0, and therefore, 
it imposes that lim,—>o = 7 . In this case, the Poisson 
equation becomes: 


y" + 2— 

r 


^2 ’ 


A5 Numerical Solutions 

The Schrodinger-Poisson (SP) system found in Eq. <§ 
does not have an analytic solution. We devote this sec¬ 
tion to finding the soliton profile for the SP system nu¬ 
merically. 


with V = log{r) a solution of the equation, for a con¬ 
stant C = -^/s. 

Thus, the SP system allows solutions with /3 = — 1, 
i.e. solutions of the form limr->o X ~ %/3/f’- Unlike the 
(5 = —2 case, these solutions will have a finite mass. 
Nevertheless, solutions with j3 = —1 are non-physical 
because they correspond to infinite total energy: 


E = f CV x)^ 4:TtA dr 
Jo 


nr 

Jo 


dr ~ 127r ( 


Finally, for /3 = 0, we must also consider the 
next terms in the expansion of x(^) small radii (i.e. 
x{r) = Uo -I- -|- ... with fdi > 0). If all the coef¬ 

ficients of terms with positive powers of r are zero, we 
recover the trivial solution x(r) = U at all radii. This 
trivial homogeneous solution corresponds to the back¬ 
ground density and the global ground state. We can 
also have non-trivial solutions for which the dominant 
term in the expansion at r —>■ 0 is a constant and the 
higher order terms have non-zero coefficients. The nu¬ 
merical solution to the boundary value problem with 
U(r —>■ 00 ) = xi’’’ 00 ) = 0 found in Appendix A5 is 

such a solution. It has limr^o — 0 > while the bound¬ 
ary conditions forbid the global ground state and en¬ 
force limr->o x” < 0 • 

The solutions with /3 = —1 and /3 = —2 are clearly 
not solutions that minimise the total energy, so they 
cannot represent the local ground state. Also, each vio¬ 
lates the non-relativistic condition, k/m <C 1 , which in 
dimensionless variables and spherical co-ordinates reads 
drX ^ 1. These solutions are not consistent limits of the 
fundamental Einstein-Klein-Gordon equations. The ex¬ 
istence of these solutions does, however, point to the 
instability of the r° flat-core (pseudo) soliton. The in¬ 
stability, however, is a relativistic one, and the SP sys¬ 
tem is not suitable to analyse it. The solutions become 
relativistic before the /? = — 1 or /3 = —2 solutions are 
ever found, and a black hole will form for finite mass 
and energy. The flat-core soliton solution is unstable to 
black hole formation when the central field value exceeds 
> O.SMpi ( Seidel fc Suen[ 1990 1. This fixes the max¬ 
imum, Chandrasekar-like, soliton mass as a function of 
rria ( |Ruffini fc Bonazzola|196^ . 

The result that the consistent, non-relativistic den¬ 
sity profiles are flat as r —>■ 0 is not surprising given 
that we were analysing solitons, which are by definition 
spatially confined, non-dispersive and non-singular so¬ 
lutions of a non-linear field theory. 


A5.1 Boundary Conditions 

We begin by noting that the SP system consists of two 
second-order ODEs with an additional unknown param¬ 
eter 7 , and thus we need to set 5 boundary conditions 
such tha t the system is uniquely determined. In Sec¬ 
tor X ~ a-s > 


tion 


A4 


^ as r —>■ 0 we found that /3 = 0 
and therefore, X^( 0 ) = 0 - In addition, we are always 
free to normalize x such that x(0) = U and then later 
restore units and use the scaling symmetry to ensure 
that x(0)^ matches the desired central density. Two ad¬ 
ditional boundary conditions can be set by imposing a 
vanishing density and potential far away from the core 
(i.e. xi’’’ 00 ) = 0 and U(r —>• 00 ) = 0). Thus, the 

first four boundary conditions arise naturally, and we 
are only left with the task of finding one more. This last 
condition arises from the fact that x^( 0 ) = x”( 0 ) = 0 
when we have a fiat core (i.e. /3 = 0). Then, by differen¬ 
tiating the Schrodinger equation with respect to r, we 
have: 

,^"' + 2^-2^ =2U'-f2(U-7)x' (A22) 


=► V\r 0) = 0 (A23) 

Thus, in order to find the soliton profile, we need 
to solve the boundary value problem: 


x" + — = 2(V - Ax 
r 

(A24) 

^ + — =X 

(A25) 

with boundary conditions: 


x(o) = 1 

(A26) 

0 

II 

0^ 

(A27) 

x(oo) = 0 

(A28) 

U'( 0 ) = 0 

(A29) 

U(oo) = 0 

(A30) 


(A31) 


Boundary value problems (BVPs) which have one 
or more of the boundary conditions set at infinity are in¬ 
herently difficult to solve numerically. One can attempt 
to replace the boundary at infinity by a boundary at 
r = 1 , but this method will be very sensitive to 

the choice of R and, in most cases, will fail to converge 
to the exact solution of the BVP. Instead, we note that 
X (r — >■ 00 ) = 0 and U (r —^ 00 ) = 0, and by assuming 
that 7 is not too small (i.e. V{R) <C I 7 I for R ^ 1), 
we can approximate the Schrodinger-Poisson system as 
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follows: 


x" + 2^ =2(V-7)x 
r 

(A32) 

x!'{R > 1) « -2 7 X (I? > 1) 

(A33) 

so at large values of R: 


x{R) oc for 7 < 0 

(A34) 

We expect 7 < 0 because x must decay at large dis¬ 
tances, whilst for positive values of 7 , x would have an 
oscillatory behaviour. If we replace this expression for 
x{R) in the Poisson equation, (A25l, we find: 

H 

(A35) 

introducing the notation U = V''. 


U'{R) + 2^^^^ oc 

R 

(A36) 

Ce-2V^R - 2 f 

(A37) 

dU dR ^ ^ , 

(A38) 

This gives the behaviour of the potential at large dis¬ 
tances: 

U{R) = V'{R) oc R-"^ 

(A39) 

V{R)<x-^ 

(A40) 

Using the approximations of Eqs. (A34l and (A40l, we 
can now replace the two boundary conditions at infinity 
(i.e. x(oo) = 0 and V{oo) = 0) by: 

xpR) = -y=^x(i?) 

(A41) 

(A42) 


for a large value of R. In practice, we will test different 
values of R and find the smallest value for which the 
solution still converges. 


A5.2 Solution of the Boundary Value Problem 


To solve the boundary value problem, convert the sys¬ 
tem of 2 second-order ODEs into a system of 4 first- 
order ODEs, by introducing the qnantities X = ^ and 
V = V'-. 


( ^ \ 

X 

V 

V V / 


(A43) 


( \ 


( X \ 

X' 


-2f+2(I/-7 )x 

V 


V 

\ y y 


1 -2^+X^ / 


(A44) 



Figure Al. Zero-node soliton profile x(r) found by solving 
the BVP system in Eq. |A25| 


with the following five boundary conditions (noting 
again that the fifth boundary condition arises due to 
the presence of the unknown parameter 7 ): 

x(o) = i, 

x'(0) = 0, 

E'(0) = 0, 

V'{R) =-V{R)/R. (A45) 

This system does not have a unique solution (y, 7 ). 
Instead, it can be viewed as an eigenvalue problem with 
unique solutions (yi, 7 i) for each eigenvalue Ai or num¬ 
ber of nodes ni. Over time, the system will relax into the 
stable state of minimum energy, which corresponds to 
a soliton with zero nodes (i.e. the ground state). Thus, 
we look not only for a solution of the BVP that con¬ 
verges, but also for the exact solution (y, 7 ) which has 
zero nodes. 

The zero node soliton profile xR) is depicted in 
Fig. |A1| We find that the zero-node solution presented 
in this figure corresponds to: 

7 = -0.692, (A46) 

V(r = 0) = -1.341. (A47) 


These values are in agreement with those reported in 
Table II of Guzman & Urena-Lopez (20041, who studied 
an equivalent system. We note that our soliton solution 
has I 7 I < 1 and |9rXl < IVr. They are consistent limits 
of the Einstein-Klein-Gordon equations. The scale pa¬ 
rameter A < 1 that takes them to astrophysical systems 
after restoring units reduces these values further and 
improves the validity of the non-relativistic, Newtonian 
limit. 


A6 Fitting the soliton density profile 

We fit the ground state soliton density profile, p = f{ar) 
for different functional forms. We begin with the form 
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Figure A2. In blue is depicted the density profile obtained 
from the zero-node soliton solution, scaled to typical dwarf 
galaxy parameters, while in red we show the fit of Eq. ( |A53[ |. 
This provides a good fit to the BVP solution up to radius 
r ^ 3.2 kpc. 


used by [Schive et al.| t |2Q14a[ > : 


f{ar) = 


1 

(1 + 


(A48) 


In Schive et al. | |2014a|b| > the soliton-NFW transition 
generally occurs for p ~ O.Olp(O). Therefore we perform 
a best fit for the value of a with r chosen so that p G 
[0.01,1]. For this range of r we find a. = 0.230. This 
form has a number of advantages. It gives an analytic 
expression for the mass integrated out to any radius. 
Assuming this form for the density, one can also find 
an analytic solution for the Newtonian potential V{r) 
from the Poission equation, and the combined solntion 
can be shown to give an accurate solution to the full SP 
system. Allowing the exponents in the fit to vary can 
improve the fit slightly at large radius, but loses these 
computational advantages. In our complete profile the 
large radius region occurs after the transition to NFW 
and improving the soliton fit here does not have any 
effect. 

We also investigate a Gaussian fit: 


f{ar) = exp 


2a2 


(A49) 


and find a = = 1.22. This form of the density 

profile also gives analytic results for the soliton mass, 
and of slightly simpler form than for the polynomial 
fit. However, we find that the Gaussian profile always 
gives a worse fit than the polynomial, and consider it 
no further. 

ISchive et al. (2014a I solved the system of ODEs in 
§ using the shooting method and the boundary con¬ 


ditions of Guzman & Urena-Lopez (20061: 


(/)(0) = a^</)(0) = 0 

(j>{r —>■ oo) = 0 
l/(oo) = 0 


(A50) 

(A51) 

(A52) 


Schive et al. 1 2014a|)’s fit for the soliton density profile 


Paoi(r) 


1.9(mB/10 ^®eV) ^(ri/a/kpc) * .3 

[l-b9.1 X 10-2(r/ri/2)2]8 ® 


(A53) 

where they define ri /2 such that Paoi(’'i/ 2 ) = Paoi(0)/2. 
This fixes the factor 9.1 x 10“^ and the value of a 
is found by normalising the central density. Restoring 
units using the scaling symmetry, as outlined in Sec¬ 
tion | 2 . 2 [ shows good agreement between this fit and 
ours. 


APPENDIX B: A NOTE ON SOLITONS 


We have considered spherically symmetric soliton solu¬ 
tions to Eqs. 0 . In this model, the non-linearity neces¬ 
sary to support solitons comes purely from gravity. Real 
fields such as axions form a class of solitons known as os- 
cillons, which, being unprotected by a charge, are tech¬ 
nically pseudo-solitons. Oscillatons are oscillons includ¬ 
ing self-gravity. The solutions we study evade Derrik’s 
theorem that otherwise forbids solitons in three spa¬ 
tial dimensions because they are time-dependent and 
include self-gravity. Scalar field solitons known as bo¬ 
son stars are formed for complex fields, where stabil¬ 
ity is guaranteed by the charge. For more discussion of 
solitons, oscillons and boson stars in various contexts, 
see e.g. [Ruffini fc Bonazzola] (|1969|; [Liddle &; Mad 


sen (19921; Liebling & Palenzuela (20121; Amin et al. 


(20121; Amin (2013); Lozanov & Amin (2014). 


We considered the ground state solitons. An iso¬ 
lated system relaxes to the ground state via the emis¬ 
sion of scalar waves to infinity, a process of ‘gravitational 
cooling’ ( [Seidel &: Suen|fl994[ [Guzman &: Urefia-Lope^ 
2006). Stable solutions are the non-linear density pro¬ 


files that form the end-point of spherical collapse in ax¬ 
ion DM halos below the threshold for black hole forma¬ 
tion. These solutions are pressure supported. By anal¬ 
ogy to stars, this justifies the use of spherical symmetry, 
as we expect the relaxation time for non-spherical per¬ 
turbations to be on the same time-scale as the relaxation 
time to the ground state. The simulations of [Schive et al.[ 
( [2014a[ ) show that these ground state solitons form and 
are stable on cosmological time scales. 

During cosmological structure formation, such 
ground state solutions will only be found locally if the 
relaxation time is shorter than the age of the universe, 
with the relaxation time being shorter for denser ob¬ 
jects. These solitons will possess a characteristic size 
related to their formation time. GDM-like structure 
will form on larger scales due to the Jeans instability 


(Khlopov et al. 19851. Furthermore, in bound states 
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16 D. J. E. Marsh and A. R. Pop 


such as halos, scalar waves may not escape to infin¬ 
ity and may contribute to the outer parts of a halo. 
Detailed questions about the non-linear dynamics can 
only be solved by simulation and will be the subject of 
future work. 
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